THE SHIFTING TECHNIQUE FOR COMPUTING THE EXTREME 
SOLUTIONS OF X + A^X-'^A = Q 
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Abstract. We propose a new way for speeding up the search of the maximal solution X+ o{X + A X~ A = Q. 
It is known that the speed of convergence of traditional approaches for solving this problem depends highly on the 
spectral radius p{X^ A). If p{X^ A) is close to one or equal to one, the iterations of traditional approaches 
converges very slowly or does not converge. Our goal is to come up with a shifting tactic to remove the singularities 
embedded in p{X^^A). Finally, an example is used to demonstrate the capacity of our method. 

1. Introduction. Consider tlic nonlinear matrix equation (NME) 

X + A'^X-^A^ Q, (1.1) 
■^'^ and Q is symmetric positive definite. We define two corresponding matrices 



where A,Q £ 

of in]) 
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Q -In 



C 



/„ 

A'^ 



(1.2) 



Note that the pencil M — XL is symplectic, i.e., it satisfies 

MJM^ = UC7 , with ./ E 



/„ 

-/„ 



and A G o{M.^ C) if and only if 1/A G <y{M., C). It is easy to see that 



M 



In 

X 



= C 



In 

X 



X-^A, 



(1.3) 



and a{X-^A) C cf{M,C)- 

For a symmetric matrix X, we use the notation X )~ Q to say that X is positive definite and the 
notation X >0 to say that X is positive semidefinite. It follows that for two symmetric matrices 
X and Y, we write X >- Y ii X -Y >- axid X >Y it X -Y >Q. 

A symmetric solution of (|l.ip is called maximal if >^ X for any symmetric solution 
X of (jl.ip . Conditions for the existence of a symmetric positive definite solution and a maximal 
symmetric positive definite solution of (|l.ip are discussed in [4 . 

Theorem 1.1. Let '4'{X) be a rational matrix-valued function defined by 



Q + XA + X-^A'^. 



(1.4) 



Then, p.ip has a symmetric positive definite solution if and only iftp{X) is regular, i.e., det'!/'(A) ^ 
for some X £ C, and ipiX) > for all \X\ — 1. 

Another necessary condition for the existence of a positive definite solution for (jl.ip can be 
described as follows. 

Theorem 1.2. // (jl.ip has a symmetric positive definite solution, then it has maximal solution 
>- 0. Moreover, is the unique solution for which p{X^^A) < 1, that is, p{X^^A) > 1 for 
any other solution X 0. Here p{ ■ ) denote the spectral radius. 

The NME arises in a large variety of disciplines in sciences and engineering. Its wide range of 
applications includes control theory, ladder networks, dynamic programming, stochastic filtering, 
and statistics (the reader is referred to [I1II3] for a list of references). Many other aspects of NMEs, 
like solvability, numerical solution, perturbation and applications, can be found in [3j H] [71 IH [lU 
[T2I [T3I [T4] and the references therein. Numerical approaches for obtaining the maximal solution 
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of (jl.ll) are mainly based on fixed-point iteration, Newton's iteration, or a structure-preserving 
doubling algorithm (SDA) [51 [71 [SJ [5] ■ However, the convergence rates of all of these methods have 
been shown to be slow while p{X^^A) « 1. Note that given a matrix norm || • |j, an algorithm 
is linearly convergent if it generates a sequence of approximate solutions {^fc} to the solution X 
such that \\Xk — X\\ < 717" for constants < cr < 1 and 7 > and quadratically convergent if 
\\Xk — ^11 < 7c^ . Our major concern is to discuss the improvement of traditional approaches 
while solving the NMEs with p{X:^'^A) sa 1. 

Our contribution in this work can be organized as follows. In Section 2, we review three 
iterative methods for finding the maximal solution of (jl.ip . In Section 3, we discuss how to shift 
eigenvalues of a general matrix pencil. A numerical example is given to demonstrate the application 
of the shifting technique. In Section 4, we summarize our results and suggest an avenue for further 
research on applications of this shifting technique for solving matrix equations. 

2. Numerical approaches. In this section we briefly review the numerical approaches for 
solving NME and the corresponding convergence rate of each method. We start our discussion 
with the fixed point iteration. 

2.1. Fixed point iteration. Let X+ be the maximal solution of (|l.ip . 
Algorithm 2.1. (Fixed point iteration for (|l.ip ) 

1. Take Xq = Q. 

2. For fc = 0, 1, . . . , compute 

Xk+i =Q- A^X-^A. (2.1) 
It has been shown in [7] that the sequence {Xk} converges to X+ and satisfies 

limsup {/\\Xk~X+\\ < piX^^Af. 

To speed up computation, Zhan [1^ incorporated the Schulz iteration [10] to provide an inversion- 
free variant method when Q = /„. This idea is then generalized in '7] for general positive definite 
Q as follows. 

Algorithm 2.2. (Inversion- free fixed point iteration for (jl.ip ) 

1. TakeXo = Q, Fo^/n/IIOIU. 

2. For fc = 0, 1, . . . , compute 

Ffe+i =rfc(2/„-Aferfc), (2.2a) 
Xk+i = Q- A^YkA. (2.2b) 

Note that Algorithm 12.21 requires more computation per iteration than Algorithm 12.11 How- 
ever, Its four matrix-matrix multiplication can be calculated in a parallel computing system effec- 
tively 0. Its numerical computation is more stable than Algorithm 12.11 due to the scheme for not 
computing the matrix inversions directly. The following result shows that the convergence rate of 
Algorithm 12. II is roughly the same as that of Algorithm 12.21 

Theorem 2.1. ^ For e > 0, the two generated sequences {Xk\ and {Yk\ of \2.2\ satisfy 

Xo>Xi>--- , Yo<Yi<-- - , 

r^+i - VII < (px-i|| + 6)2||r, - x-^ll, 

and 

\\x,-x+\\<\\Arm-x-'\\. 

It follows that the convergence rate of iteration (|2.1|) or ()2.2p is highly related to the spectral 
radius p{X^^A). Also, we see that Algorithm l2 . 1 1 and 12 .21 are linearly convergent when p{X^^A) < 
1 and the convergences slow down when p{X^^A) is close to 1. In these situation, an alternative 
approach, Newton's method, is recommended. 
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Then the Frechet derivative TZ^ : 5" 



2.2. Newton's method. Let be the set of positive definite matrices in R"^" and 5" 
be the set of all symmetric n x n real matrices. Corresponding to (|l.ip . we define an operator 
n-.V" 5" satisfying 

n{X) = -X + Q - A^X-'^A, X-^O. (2.3) 
5" of 7?. at X is given by 

n'xiZ) = -Z + A'^X-^ZX-^A, ZyO. (2.4) 
By (|2.3p and (|2.4p . we obtain the Newton step for the solution of (jl.ip is 

= Xk-i - {nx,_y'n{Xk-i), fc - 1, 2, . . . . 

Combining (|2.3p and (|2.4I) . we have the algorithm of Newton's method for p.ip . 
Algorithm 2.3. (Newton's method for (|l.ip ) 
J. Take Xq ^ Q. 

2. For fc = 0, 1, . . . , compute Lk ~ X^^^A and solve 

— LjXkLk = Q — 2LIA. 

It is known that in each iteration, the computational work for Newton's method is about 
15 times larger than that for the fixed point iteration. In order to have a better picture of the 
advantage of applying Newton's method, we include the convergence result discussed in [7] as 
follows. 

Theorem 2.2. // (jl.ip has a positive definite solution, then Algorithm \2.3\ determines a 
nondecreasing sequence of symmetric matrices {^fc} for which p{Lk) < 1 and lini Xk — X+. 

k^oo 

Moreover, if p(X^^A) < 1, the convergence is quadratic and if p{X^^A) — 1 and all eigenvalues 
of X^^A on the unit circle are semisimple, the convergence is either quadratic or linear with rate 
1/2. 

Theorem 12 . 21 states that if p{X_^_^A) = 1, convergence is guaranteed under a addition assump- 
tion, i.e., all eigenvalues of X^^A on the unit circle are semisimple. In [8], Lin and Xu investigate 
the approach of applying the structure-preserving algorithm (SDA) for solving (jl.ip without any 
assumption on the unimodular eigenvalues of X^^A. 

2.3. SDA. The fundamental idea of SDA is the doubling transformation. In this sense, we 
begin with the discussion of the doubling transformation. 

Let Ai — XC be a matrix pencil consisting of two matrices 



M = 



A 

Q -In 



and C = 



-P In 

A^ 



(2.5) 



with Q, P ^ 0. This is the so-called the second standard symplectic form (SSF-2). For any 

C 

-M 



solution [A^*,£*] within the left null space of 
transformation 



, define M := Mi,M and C := The 



M-XC^M-XC 
is called a doubling transformation. Assume further that 



A{Q-Py^ 

-AT(g-p)-i /„ 



and Ci, 



In 




-A{Q 



P) 
P) 



We see that by direction computation 
M:=M*M = 



A 

Q -In 



and L := = 



-P In 

A^ 
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where 



A:=A{Q-P)-^A, 



A'^iQ - P)-^A and P P + A{Q - Py^A'^. 



(2.6) 



This impHes that if Q-P ^ and Q- ^^(Q-P)-^ >z 0, then {M, C) is again a SSF-2 form [S]. 
Based on formulae (|2.6p . we then have the foUowing algorithm, SDA. 
Algorithm 2.4. (SDA for pTT]) ) 

L Take Ao = A,Qq^ Q, Pq = 0. 

2. For fc = 0, 1, . . . , compute 



Ak+i = Afe(Qfc - Pky^Ak] 
Qk+1 = Qk - AjiQk - Pk)-^Ak 
Pk+i ^Pk + AkiQk - Pkr'Aj. 



(2.7a) 
(2.7b) 
(2.7c) 

is well-defined, that is. 



Below we quote from [S] Theorem 2.1] to guarantee that Algorithm [ 
the difference — >- 0, for all k. 

Theorem 2.3. JEI Let X > he a solution of (dH]). Define S ^ X-'^A. Then the sequences 
{Afc, Qk, Pfe} generated by Algorithm \2.4\ satisfy 

1. Ak = {X ^ Pk)S^' ; 

2. 0< Pfe < Pfc+i < X and Qk - Pk = {X - Pk) + Al{X - Pk)-^ Ak > 0; 

3. X< Qk+i <Qk<Q and Qk-X = {S^f{X - Pk)S^" < {S'^fxS^\ 
Moreover, we have 

1- \\Ak\\2<\\Xh\\S^'h; 

2. \\Qk-Xh<\\XU\S^''\\l 
From Theorem 12.31 we know that if p{X ^A) < 1, then the SDA is quadratically convergent. 
If p{X-^A) = 1, Chiang et al. in g] proved that Algorithm O for NME is linearly convergent 
with rate 1/2, without any assumption on the unimodular eigenvalues of p{X~^A). In this case, 
we are interested in exploring a strategy of shifting unimodular eigenvalues of X~^A so that the 
convergence of above algorithms can be speeded up. 

3. The Shifting technique. Assume that AT )^ is a solution of (|l.ip . Then the solution 
X is highly related to the generalized eigenspace of the pencil 



M-XC 





-Ir, 



X 











That is, if X ;^ is a solution of (|l.ip if and only if X satisfies (|1.3p . 

Corresponding to (|1.3I) . let us focus on the discussion with the shifts of eigenvalues of the 
matrix pencil Ai — XL. To begin with, we consider a single shift of an eigenvalue of the matrix 
pencil M - XL. 

Lemma 3.1. Let Ai — XL be a matrix pencil with Mv = XqLv for some nonzero vector v. If 
r is a vector with r'^v = 1, then for any scalar Xi, the eigenvalues of the matrix pair 

M-XL = M + iXi- Xa)Lvr'^ - XL, 

consist of those of Ai — XL, except that one eigenvalue Xq of Ai — XL is replaced by Xi. 
Proof. Since Aiv — XqLv and r'^v = 1, we have 



Note that (M - XC)v ^ XqCv - XCv = 

det{M - XC) = dct{M 
= det{M 
_ Al - A 
" An - A 



Aiv = XiLv. 

- (Xq — X)Cv. Also, for any A 7^ Aq, we see that 

- XC) det(/„ + (Al - XQ)Lvr'^{M - XC)~^) 

- A£)(l + (Al - Ao)r^(X - XCy^Cv) 

det(X - XC). 



(3.1) 



(3.2) 
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Thus, the theorem foUows from p.ip and p.2p . 
□ 

Similar to the proof given above, we then have the following result that k eigenvalues of Ai — XC 
are shifted simultaneously. 

Theorem 3.1. Let Ai — XC be a matrix pencil with eigenvalues Ai, • • • ,Xk that satisfy 

Mvi = XiCvi,- ■ ■ ,Mvk = XkCvk 

for some nonzero vectors wi, • • • , Wfe G C". Suppose V — [wi, • • • ,Vk] G C"^*^, A = diag{Xi, ■ ■ ■ ,Xk) 
and A = diag{Xi, • • • , Afc) G C*^^*^ for some Ai, • • • , Afc G C. If Ri and R2 are two matrices in 
C"^'= such that 

RJV = A- a and RJV = 0, (3.3) 

then the eigenvalues of the matrix 

M-XC:= {M + CVRJ) - X{£ + MVRJ), 

consist of those of M~XL, except that eigenvalues Ai, • • • , A^ ofM — XC are replaced by Xi, - • • , Afc. 
Moreover, Ai,-- - , Afc are eigenvalues of M — XC corresponding to the eigenvectors vi,--- ,Vk, 
respectively. 

Proof Since MV = CVA, we have {M - X£)V = CV{A - A/„). Thus, {M - XC)-^CV = 
V{A — A/„)^^. This implies that for any A 7^ Ai, • • • , Afc, the determinant of — XC is 

det{M-XC) = dct{M - XC) det(/„ + {CVRJ - XCVAR^){M - XCy^) 

= det{M - XC) det(/fc + rJ{M - XC)-'^CV - XARJ{M - XC)-'^CV) 

1 1 1 



det(X - XC) dot ( 4 + (A - A)diag 
det(X - XC) 



Ai — A ' A2 — A ' Afc — A 
Ai — A\/A2 — A\ /Afc — A\ 



Ai — Ay\A2 — Ay y Afc — Ay 
Also, by p.3p . we have 

MV = CVA = CVA. 

This completes the proof. □ 

Note that the matrix Ri can be obtained by using the Gram-Schmidt process to the column 
vectors of V and R2 can be obtained from the vectors in the orthogonal space of the space spanned 
by the column vectors of V. We use the following example to demonstrate an application of the 
shifting technique discussed above. 

Example 3.1. Assume n = 1. Then (jl.ip can be written as 

X + — = (?, (3.4) 

X 

where a, g G M and g > and the corresponding matrix pencil is denoted by 



Ml - XCi 



a 

q -1 



1 
a 



(3.5) 



Let A = — 4a^ be the discriminant of p.4p and T{X) = det {Mi — XCi) = —aX^ + qX — a be the 
determinant of (|3.5|) . Corresponding to p.4p . we define the function il^i{X) such that 



(3.6) 
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It is clear that •ipi{X) is regular. Let A = e'^ for every 9 ^ R. Substituting this A into 
we obtain ipi{e^^) = 2acos{9) + q. This implies that ■0i(e*^) > for every 9 E M. if and only if 
A > 0. Upon using Theorem ] 1. 11 we know that (j3.4l) has no symmetric positive definite solution, 
provided A < and has a maximal solution > 0, provided A > 0. In particular, a;+ > \a\, since 
pix^^a) < 1 (see Theorem \1.2\) . 

Now we are ready to rewrite Algorithm (|2.4|) corresponding to p.4p as follows. 

Algorithm 3.1. (SDA for solving (|X4|) ) 

1. Take ao = a, qq = q, po = 0. 

2. For fc = 0, 1, . . . , compute 



at 



<7fc+i 



qk -Pk 

qk — afc+i; 



pk+1 = Pk + flfc+l- 

// the maximal solution x+ > Q of p.4p exists, by Theorem \2.3\ we have 

Qk = {x+ -Pk){ — Y , 



(3.7a) 

(3.7b) 
(3.7c) 



Qk-Pk^ {x+ -pk) + 



qk-x+ = {x+ -pk){ — f 



Pk 



Observe further that if X-^- — \a\, we obtain q = + — 2|a| > by substituting X-^. into p. 4 
It follows that A = and any positive solution satisfies 



= 0. 



Without loss of generality, assume a > 0. By induction and p.7|) . it is easy to see that 

(2'= + l)a (2'=-l)a 



a 
2^' 



qk 



2k 
1 



qk 



-, for k = 0,l,. 



limsup {/\qk - a\ 



This means that the sequence qk converges linearly to a with rate i. On the other hand, T[X) 



-a(A — 1)^ implies that p.Sp has the the right eigenspace Ei = Ker{Mi — lCi\ = Span{ 



} 



corresponding to the eigenvalue 1. Let vj — [l,a]. By Theorem \ 3.1[ define two corresponding 
vectors rj — [r — 1,0], for some constant < r < 1, and rj = [0,0] such that the new matrices 
A4i and Ci is given by 



Ml Ml + CiVirJ 



£i := Ci +MiVirJ 



ar 

a{l + r) -1 

1 
a 



By direct computation, we see that the eigenvalues of the new matrix pencil Mi — XCi are r and 
1 and the eigenvectors can be chosen as [l,a]^ and [l,ar]^ , respectively. Let vj = [l,ar]. Using 
the same constant r given above, we define two vectors rJ = [■^ — 1,0] and fj = [0,0] such that 



M2 ■= Ml + Civifj = 
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£2 := A + Mivifj = 

it is worth noting that structure the pencil M2 — A£2 is indeed SSF-2. Also, the eigenvalues of the 
matrix pencil M.2 — ^£2 are r and - with eigenvectors [1, and [1, ar\^ . Note that the maximal 
solution xj^ of 

a;+ — = ar + -), (3.8) 

X r 

is x+ — ^ and ^ = r < 1 . This implies that if we apply Algortihm \8.1\ to find the maximal 
solution of (j3.8p the convergence rate of the iterations is quadratic. By the continuity dependence 
of the solution of the NME, the maximal solution of p.4|) can be obtained by taking r — )■ 1~ so 
that X-\- approaches x+, the maximal solution of p.4|) . 

4. Conclusion. In this work, we provide an approache of shifting eigenvalues of general ma- 
trix equations. Our goal is to remove the singularities happened while solving NME. Currently, we 
have not applied this method to a much more general NME or any other nonlinear matrix equation. 
We believe this research would propose an avenue for speeding up the numerical approaches for 
solving nonlinear matrix equations. 

Acknowledgment. This research work is partially supported by the National Science Council 
and the National Center for Theoretical Sciences in Taiwan. 



1 
a 



REFERENCES 

[1] W. N. Anderson, Jr., T. D. Morley, and G. E. Trapp. Positive solutions to X = A- EX'^B* . Linear Algebra 
Appl, 134:53-62, 1990. 

[2] Chun-Yueh Chiang, Eric King-Wah Chu, Chun-Hua Guo, Tsung-Ming Huang, Wen- Wei Lin, and Shu-Fang 

Xu. Convergence analysis of the doubling algorithm for several nonlinear matrix equations in the critical 

case. SIAM J. Matrix Anal. Appl, 31(2):227-247, 2009. 
[3] J. C. Engwerda. On the existence of a positive definite solution of the matrix equation X -f A^X~^ A = I. 

Lin. Alg. Appl, 194:91-108, 1993. 
[4] Jacob C. Engwerda, Andre C. M. Ran, and Arie L. Rijkeboer. Necessary and sufficient conditions for the 

existence of a positive definite solution of the matrix equation X + A*X~^A = Q. Linear Algebra Appl., 

186:255-275, 1993. 

[5] G. H. Golub and C. F. Van Loan. Matrix Computations, 3rd ed. The Johns Hopkins University Press, 1996. 
[6] C.-H. Guo. Convergence rate of an iterative method for a nonlinear matrix equation. SIAM J. Matrix Anal. 
Appl, 23, 2001. 

[7] C.-H. Guo and P. Lancaster. Iterative solution of two matrix equations. Math. Comp., 68:1589—1603, 1999. 
[8] Wen- Wei Lin and Shu-Fang Xu. Convergence analysis of structure-preserving doubling algorithms for Riccati- 

type matrix equations. SIAM J. Matrix Anal Appl, 28{l):26-39, 2006. 
[9] B. Meini. Efficient computation of the extreme solutions of X -|- A*X~^A = Q and X — A*X~^A = Q. Math. 

Comp., 71:1189-1204, 2002. 

[10] Gnther Schulz. Iterative berechung der reziproken matrix. ZAMM - Journal of Applied Mathematics and 

Mechanics / Zeitschrift fr Angewandte Mathematik und Mechanik, 13{1):57— 59, 1933. 
[11] J.-G. Sun and S.-F. Xu. Perturbation analysis of the maximal solution of the matrix equation X + A^ X^^ A = 

P, II. Lin. Alg. Appl, 362:211-228, 2003. 
[12] S.F. Xu. Numerical methods for the maximal solution of the matrix equation X + A^X~^A = /. Acta 

Scientiarum Naturalium Universitatis Pekinensis, 36:29-38, 2000. 
[13] X. Zhan. Computing the extremal positive definite soluions of a matrix equation. SIAM J. Sci. Comput., 

17:1167-1174, 1996. 

[14] Z Zhan and J. Xie. On the matrix equation X + A^X'^A = I. Lin. Alg. Appl, 247:337-345, 1996. 



